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ABSTRACT 


Evaluating natural frequencies and critical speeds of 
rotors is an important problem since these are integral 
components of many machines. 

Transfer Finite Element Method has been used in the 
present work for determining natural frequencies and critical 
speeds of rotors taking into account shear deformation, 
rotary inertia, and gyroscopic action effects. This method 
offers the advantage of reduced computer memory requirements 
vis-a-vis standard Finite Element Method. A variety of 
situations are considered and results are found to be 
satisfactory. 



CHAPTER 1 


INTRODUCTION 


Dynamics of rotating shafts, rotor dynamics for short, 
is an important area in the field of mechanical engineering 
by virtue of the fact that rotors are an integral part of all a 
turbo-machines and also find use in many other machines. 

Modern machines frequently use ’high-speed’, also called 
’flexible’, rotors. These are rotors running at speeds above 
critical which marks the phenomenon of unacceptably large 
transverse vibrations and whirling, A complete understanding 
of the dynamics of such a rotor would involve an understanding | 
of the shaft behaviour in terms of critical speeds, behaviour ii 

ISt ... i 

post-critical speed regions, the stresses developed, interactioi 
between supports and the shaft, the effects of material | 

damping and resistance between shaft and bearings, etc. A I 

complete analysis of these phenomena is obviously beyond the | 
scope of a work such as the present one,. 

Since the 1870’s when flexible shaft theory first made ; 
its appearance in connection with the then novel turbine I 

techniques [1] much work has been done on various aspects of 
rotor dynamics. Research on the phenomenon of shaft vibrations 
was initiated by Rankine and Laval in the last quarter of the 
nineteenth century [1] . A.Stodola, an eminent turbine expert,; 
carried out investigations on many fundamental phenomena in the 
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early part of this century. More recently, Dimentberg [ 1 ] 
and Tondl [2 ] have published excellent works on rotor dynamics. 
Rao [3] has discussed, amongst other things, shaft and fluid- 
film bearing interaction at length. 

The present work focusses attention on a specific area - 
viz. finding the natural frequencies and critical speeds of 
rotors under various conditions. The specific problems address 
are: 

i) uniform rotors with different boundary conditions 

ii) non-uniform rotors 

iii) rotors with in-span supports 

Problems of this nature have, of course, been solved 
earlier using different methods. (Section 1.1 outlines 
some of the more important ones). The present work applies 
the Transfer Finite Element technique [4] for solving the 
problems mentioned above. Some notes on the Transfer Finite 
Element Method are given in Section 1.2. Different FORTRAN 
programs have been developed to obtain values of critical 
speeds for different rotor systems. 

It is highly advantageous to work with nondimension- 
alised variables of space and time for it makes the results 
very general. The results become applicable to shafts which 
are different physically but have the same nondime ns ional 
parameters. The governing differential equations have therefoj 
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nondime ns i on al quantities. 

Methods of Analysis 

The methods generally used for finding natural frequen- 
cies of transverse vibrations can be grouped broadly under the 
following heads: 

i) Exact (closed-form) solutions 
ii) Approximate methods. 

Some of the popular methods in these categories are 
mentioned here, 

i) Exact Solutions: The earliest efforts at determinii 
natural frequencies of rotors treated the shaft as an Euler 
beam [5 ] . Indeed, if the effects of rotary inertia, shear 
deformation, and gyroscopic action of discs and rotor elements 
are neglected, the problem does reduce to one of an Euler 
beam. Making these assumptions, then, it is possible to obtain 
an exact expression for natural frequencies from the governing 
differential equations of the system. 

Dimentberg [1] has obtained closed form solution for 
a simple-simple rotor considering all the effects mentioned 
above. There are, however, two major limitations: (i) result 
can be obtained only for simple rotor configurational geometri 
and (ii) determining of the frequency equation is a very 
tedious job specially if the end conditions are not simple. 
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Gyroscopic action of a disc mounted on a shaft is 
considered in many texts [1,3, 6, 7] . In all these cases, 
however, the rotor itself is considerad massless and effects 
of rotor elements are ignored. 

Patnaik and Mallik [ have investigated the natural 
frequencies of a rotor with lumped inertias considering the 
mass of the rotor also. They consider a rotor on multiple 
supports. However, the shaft and bearing structure is assumed 
to be periodic which limits the gaierality of their analysis. 

The scope of exact solutions is, therefore, severely 
limited. 

ii)- Approximate Methods: Included under this are a 

number of techniques. Applied to problems in vibrations, some 
of the popular ones are: 

a) Rayleigh-Ritz method 

b) Galerkin method 

c) Transfer Matrix method 

d) Finite Element method 

Rayleigh-Ritz and Galerkin methods [13] require one 
to assume a series solution for the system. All basis 
functions in this series solution should satisfy the 
boundary conditions of the problem. The difficulty lies in 
choosing these basis functions to start with. As pointed 
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out in [7] the method calls for a certain degree of skill 
on the part of the user in its application. 

Transfer Matrix Method (TMvl); Pestel and Leckie [9] 
offer a first-rate treatment of the subject. Pilkey and 
Chang [10] have also applied the method to a wide variety 
of problems. Their text is in the nature of a handbook and 
contains exhaustive tabulation of transfer matrices for 
different systems. 

In applying TIvIM one first needs the solution of the 
governing differential equation [9,10]. From this, one 
constructs the transfer matrix to relate state vectors on 
neighbouring stations; The state vector is transferred from 
one station to the next by ’field transfer matrices', brrerfly 
'field matrices'. For points across which one or more state 
variables are discontinuous (like thin discs, in-span supports) 
state vector is transferred across this point by a 'point 
transfer matrix', briefly 'point matrix'. The state vector 
can thus be transferred across the entire rotor thereby 
relating state vectors at the rotor boundaries. Boundary 
conditions of the problem are then imposed and natural 
frequencies determined. 

Although this is a powerful method, the limitation 
lies in the fact that many times one cannot find a general 
solution to the differential equation. Moreover, even when 
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this is possible, the expressions can be quite cumbersome 
and derivation of the transfer matrices a tedious task. 

Finite Element Method (FEM): This is a powerful technique 

described in various texts [12, 20] , 

FEM does not suffer from the limitation pointed out 
for above. Here generally a polynomial solution over 
the finite element is assumed, A matrix form of equations 
between state variables, valid for the finite element, is 
derived . Matrix aquations for various elements are next 
assembled together. 

This process of assembly leads to an increase in the 
order of matrices being handled. For a complicated system 
having many elements these matrices can become very large 
leading to large computer memory requirements. This drawback 
apart, FEM is a very powerful technique and is being widely 
used to solve a host of complex problems in a variety of 
fields in engineering. 

1.2 Notes on the ‘Bransfer Finite Element Method (TFEM) ; 

Gupta [4] has developed this elegant technique 
utilising ideas from both TMM and FEM. The method ccanbines 
the advantages of both - reduced demand for computer memory 
without losing the polynomial solution advantage of FEM, 
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The matrix equations of the finite element obtained 
in the usual FEM are rearranged to yield transfer matrix. 

This relates state vectors at the neighbo^^uring nodes as in 

im. 

The TMA technique of multiplying these transfer matrices 
is used to obtain the overall transfer matrix [ OT] between 
state vectors at the boundaries of the system. A determinant 
of residue is also obtained on similar lines by imposing 
boundary conditions. Hence ’assembly’ for the entire system 
is different from the standard FEM and results in saving on 
computer memory, 

Subramanya [11] has used this technique for the 
specialised topic of flexibility analysis of pipes. 



CHAPTER ■ 2 


FORMULATION OF GOVERNING EQUATIONS 


Earliest efforts at analysis of continuous systems of 
rotors made simplifying assumptions of neglecting rotary 
inertia, shear deformation, and gyroscopic effects. These 
simplifications lead to the classical Euler-Bernoulli formula- 
tions, However, results obtained from this formulation do not 
match known experimental values of natural frequencies and 
critical speeds, specially for higher modes and for short, 
stubby rotors. Lord Rayleigh and S, Timoshenko went into 
deeper analysis and showed that the effects of rotary inertia 
and shear deformations can be very significant in these 
cases. For rotors with large diameters and high rotational 
speeds, gyroscopic action of elements can also be significant. 

Rotary inertia effect arises due to the fact that 
rotor elements not only move laterally in the plane of 
vibration but also rotate about an axis normal to this plane 
when the rotor bends. Again, shear deformation effect arises 
from the fact that shear forces on element faces cause defor- 
mation of the element and affect total slope. Total slope 
now is due to both bending moment and shear force. Both 
rotary inertia and shear are known [21] to decrease natural 
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frequencies vis-a-vis values obtained from classical Euler- 
Bernoulli theory. 

Gyroscopic action occurs in rotors spinning about an 
axis passing through their geometric centre. Because of this 
rotation there is an angular momentum vector associated with 
each element. Transverse movement of the rotor causes a change 
in the angular momentum resulting in a gyroscopic couple being 
produced. As Rao j^S] shows, this couple tends to straighten 
the shaft, i.e, it tends to realign the angular momentum vector 
with the centre line of the rotor. It therefore makes the 
rotor stiffer and increases values of natural frequency. 

Critical speeds of a rotor are rotational speeds equal 
to the natural frequencies of transverse vibrations of the 
system. These are so called because all rotors, no matter how 
well balanced, have some amount of residual lirbalance and this 
causes resonance when the rotor spins at a speed equal to its 
natural frequency. Whirling can become very large and even 
cause damage in such situations; 

Equations are presented here for rotors taking into 
account rotary inertia, shear, and gyroscopic effects. Two 
approaches are presented — (i) two governing equations and, 
(ii) these combined into a single governing equation for 
the special case of rotor-segments with uniform cross-sectional 
properties. Some comments on these two approaches are presented 
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later in this chapter; 

Equations are presented first in the dimensional form 
and then non-dimensionalised using characteristic values of 
length and time. 


2,1 Basic Equations for Rotor : 

These have been derived by various authors [1f2] and 
are being presented here for sake of completeness. 


Rotor motion does not occur in one plane. Instead, 
the rotor spins about its axis as well as processes about 
the bearing axis. It is thus necessary to consider the motion 
in two mutually perpendicular planes passing through the bearing 
axis. Coordinate system used is shown in Fig, 2, 1(a), Consider 
a thin section of length A s on the rotor (referred to as ’ eleme 
in this chapter but not to be confused with finite elements) 
in planes s-x and s-y. 


Figure 2.1(c) shows the element in plane s-x. 
Figure 2.1(d) shows the details of slope of the element 
■W y is slope due to bending moment only while 3 x/3 s is 
slope due to bending moment and shear. 


total 


From solid mechanics [13] we have the relations: 



El 

k' GA ( 


d'S 


3 s 
3 s 



( 2.1 . 
( 2 . 1 , 



X 

A 


] >s 


X 



(a) Coordinate System 


X 

A 


T 

X 

u 


"y/ 

X S, 


As 

N ^ 8s 

o 


’) 


(b) Angles subtended by angular 
Momentum Vector 

centre-ltne 


Mv'*' ^My- As 
^ 8s 


->s 



(c) Rotor section in plane s-x 


(d) Details of slope of rotor 
section in plane s-x 



(e) Rotor section in plane s-y {f ) Details of slope of rotor 

section in plane s-y 


4 CCrTlHMQ 



12 


where k* is the appropriate form factor and is 0.9 for circular 
cross-section. 


From force balance in x-direction 
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First term on RHS of above equation is rotary inertia 
term while the second one is contribution from gyroscopic 
action of the elemental disc. 

Similarly, for plane s-y, we have Fig. 2. 1(e), (f). 
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®“x 2 2 

= mk^ ^ + 2mk"^0)-^ (2.1.2( 

d s y ax 


2*2 Two Equation Approach ; 


Straightforward substitution of relations for M and S 
into the last two of equations (2,1,1) and (2,1,2) yield f'or 
plane s-x 
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Multiplying equations (2.2,2) by i = -1 and defining 

canplex displacement as z = x + iy, complex slope as W =W y -iW ^ 
one can combine the above equations to get: 

mi- - -I; [k'GA(z' - 5 ) ] = 0 (2_2.3) 

-4- (El 5^') + k’GA(zW ) - mk^f’ + 2imk^ = 0 
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where primes denote derivatives w.r.t. s and dot w.r.t. 
time, t. 

Equations (2.2.3) are the coupled governing differential 
equations for a rotor. 


In order to obtain non-dimensional form of the above 
equation, for the present work the special case of constant 
cross-sectional properties is considered. The two equation 
approach has been used only for comparison of results with 
one-equation approach for this specific case. As pointed out 
earlier, other problems have been solved using the single 
governing equation approach. 

For constant cross-sectional properties, non-dime ns ion- 
alised equations are (See Appendix A) 




(Z -5- ) - 


Mc.Re^.Z = 0 


9 ^' + — J — _ (z»«) - Re^ (?* - 2iC f) = 0 (2.2, 


3S 


Me .Re 


2.3 Single Equation Approach : 

Equations (2,2.3) above can be combined to yield a 
single fourth order governing differential equation for the 
special case of constant cross-sectional properties . 
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Assuming that k*GA is a constant, we obtain from 
relations ( 2 . 1 . 1 b,c) the expressions 

m 9 X 

m 3 

k‘GA 3523^2 

and from ( 2 . 1 . 2 b,c) the expressions; 
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Assuming again that El is constant, differentiating 
equations ( 2.1 .1 d) and ( 2 . 1 . 2 d) once w.r.t. s and using 
expressions ( 2 . 3 . 1 ) and ( 2 . 3 . 2 ) we obtain equations for the 
two planes. These can be combined as for the two-equation 
approach to yield the single governing equation in z; 
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Equation (2.3,3) (as also equations (2.2.4)) takes 
into account rotary inertia, shear effects, and gyroscopic 
effects. 

Timoshenko beam and Euler beam equations are special 
cases of equation (2.3.3). The former requires us to neglect 
gyroscopic effects and the equation can be established by 
taking Q = 0 in this equation, Euler beam formulation neglects 
all the above mentioned effects and hence the equation can be 
obtained from equation (2.3,3) by retaining only the and z* 
terms. 

In the non-dimensional form equation (2,3.3) is 
(Appendix A): 


TV o •«ti A • •• 

Z - Re (1+Mc) Z + 2iA.Re Z + Me. Re . Z + Z 


-2i52.Mc.Re^ Z = 0 


(2.3.4) 



Consider. the rotor shown in fig, 2, 2 below: 
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The three segments have cross-sectional areas 
and A^ with corresponding diameters d^ , d^ and d^. E, G and 
k’ are considered same for all sections. 


Equation (2,3.3) is applicable to a shaft section with 
uniform cross sectional properties. This equation is therefore 
applicable to sections 1,2 and 3 of the shaft with suitable 
values of cross-sectional variables. Thus for jth section 
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To non-dimensionalise this equation, one can use overall 
shaft length L, and characteristic time T^ = ( mL'^/EI) based 
on characteristic values of m and I -of any section. Denoting 
these characteristic values by non-subscripted quantities, 
one has from equation (2,4.1): 
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Denoting d./d by Rd. (relative diameter of section j 
with respect to characteristic diameter, d), we have; 

o TV oof * * " i.'” ■ 

Rdj. - Re .Rdj (1 + Me) Z - 2i: 0 Z 
' + *Z* + Mc.Re^.Rd? * Z*- 2i D . Mc.Re'^.Rd? *Z = 0 

%J %J 

(2.4.2] 


where Re = k/L is based on characteristic k. 

The above equation is, for the present work, the most 
general form of the governing differential equation for trans- 
verse vibrations of a rotor. Putting Rd^ = 1.0 for all sections 
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would yield, the equation for a uniform rotor (ref. equation 
2 . 3 . 4 ). 



CHAPTER 3 


FINITE ELEMENT FORMULATIONS 


In the course of the present work, it was observed that 
using two coupled differential equations for obtaining natural 
frequencies and critical speeds leads to problems of poor conver- 
gence. This happens because equations become ill-conditioned. 
Reduced numerical integration using Gaussian quadrature has to 
be carried out to overcome this problem [19,20 ] . It may 
mentioned that if one is interested in situations where cross- 
sectional area of the shaft varies continuously along the 
shaft length, one shall necessarily have to use a two-equation 
approach for it does not seem possible to combine the two 
equations to yield a single governing differential equation 
for such a case. However, if instead one is interested in 
situations where shaft section changes are discrete with cross- 
sectional area being constant for a given segment, one can 
use a single governing equation approach. 

The two equation approach is also the more useful one 
when one is faced with situations where both ends of the beam/ 
rotor are built-in (or fixed). This is so because by using 
this approach one can satisfy both the geometrical boundary 
conditions exactly. 
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Details of the finite element formulation have been 
provided here with reference to the single governing equation 
of fourth order. Formulation for the two-equation approach 
is presented in Appendix B. 

From the previous chapter it is seen that for a rotor 
equation (2,4,2) is more general and differs from equation 
(2,3.4) in only some constant factors to account for differences 
in cross-sectional properties. Equation (2.4.2) is therefore 
used as the basic one for deriving finite element fomulations 
here. This equation is presented again, with minor modifica- 
tions, for ease of reference: 

C^Z - [ (1 + Me) Z - 2 i Q Z J + Z 

+ C 3 . *Z * - 2 i Q . C 3 . *Z * = 0 (3.1 ) 

where 

C-i = Rdj, = (Re * C 3 = Mc.Re'^.Rd? . 

Consider the rotor divided into number of finite 
elements, as in figure 3.1 
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Node No,: 1 2 3 j (j+D (j+i) 


n Co+0 



Fig.3,1:(a) Rotor divided into finite elements 
( b) Typical finite element 


A typical element is also shown in the figure. Local 
coordinate along the length is g and element length is h. 

In non-dimens ionalised form, we have T) = $ /L and H = h/L. 

The governing equation applicable to the element is 

then: 
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where 

Z = Z (T], T) . 

Assume a polynomial solution over the element, 

(ti,T) - [n^ N^ ... N^J[z^ Z2 *** ^rj^ 

= LnKz}^"®' (3.3) 
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where are shape functions dependent only on and are 
variables at nodes of element and are functions of time only. 
Since this is an assumed solution, it will not satisfy the 
governing equation exactly. Substitution of expression (3.3) 
in equation (3.2) leaves a residue over the element. 
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(3.4) 


This residue needs to be minimised. Different techniques 
are available for this. Here the most successful Galerkin 
method [i6] has been used. Thus one gets 

^ (e) 

f N. . Res^®^ . dY) = 0 

i~1.2......r ( 3 . 5 ) 

Equation (3.5) is integrated by parts, till the derivatives 
are reduced to the lowest order possible. Using expression 
(3.3) this becomes; 
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No comment has so far been made on the number of terms 
in the polynomial approximation (3.3) and on the nature of shape 
functions N^, Various finite element models using different 
shape functions have been used in the literature for the dynamic 
analysis of Timoshenko beams, Thomas and Abbas [14(] have given 
a comparison of different models used by various investigators 
and have also presented a better finite element model. Earlier 
models, they have claimed, did not represent the Timoshenko 
beam accurately because all boundary conditions could not be 
applied satisfactorily. Their model is an element with ei^t 
degrees of freedom requi ring eight nodal coordinates , They have 
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assumed cubic polynomial approximations for total deflection 
and bending slope of the element. No justification has been 
given for this choice of higher order element. 

Gupta [4] and Singh and Singh [15] have used a rational 
approach of Huebner [16] for deciding the number of terms in 
the interpolating polynomial. In both cases finite element 
equations are systematically and logically derived from the 
governing differential equations of the system using the Galerkin 
method and the order of the approximating polynomial is obtained 
from the conditions of compatibility and completeness as outlined 
by Huebner [16] . This approach is followed here. 

The finite element equations (3.6) have been derived 
using the Galerkin method from the governing differential 
equation (2.4.2). The highest order of the derivative appearing 
under the integral sign in equation (3.6) is two which requires 
a compatibility of Z and first derivative. The highest order 
derivative in the same equations is three and hence completeness 
of Z and derivatives upto third order is required. Therefore, 
the approximating polynomial needs to be of third order. 

Thus the finite element model for this problem can be a two 
noded element with Z and Z* as degrees of freedom at each node. 
Shape functions satisfying the above requirements of compatibility 
and completeness are [12] : 
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N = [ 1 


+ »£ 
h 2 


22 : 

H 


X]- + 


It 


3ti^ 2t£ 




H" 


i£ + ni -j 

H ^ 2 ] 
(3.7) 


and nodal variables are: 


= [Z^ 


Z 


1 


1 7 

^2 3 


(3.8) 


Using expression (3.7), various matrices of equation( 3.6) 
can be evaluated and this equation written as 


[K'] + C 2 (1 + Me) 

; (ne) .. (ne) 

- C 2 2i Q[MB] (Z) +[M] {ZJ 

.... (ne) ... ( ne) 

+ C 3 [M] tZ } - 2i Q C 3 [M ]{Z } 




< 


(C, 


- c. 


3.(e) 


d^^Z 


dn 


— — C 2 ( ( 1+Mc) 


3-(e) 


- 2 i 


, d^Z 


2 ,(e) 


2y) 




)) 


T 1=0 


d¥®3 

^ 7 ^ 


j3_(e) 

(C d Z_ _ c (d+Mc) -2i 

^ dT)'^ dr)dT 


in= 0 


i^ 7( ® ^ 3 ^ ^ 


”!n 


2,(e) 


r ^ 
d-^ 


>n= H 


( 3 . 9 ) 
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where 

[K] = f (n"} Ln"J dti 

o 


[M] = I (n}1_n J dn 

O 

H 

[MB] = j {N’llN'J dt) 
o 


are square matrices of order 4x4, The elements of these matrices 
with LnJ as in expression (3,7) are given in Appendix C, 

We can expect the system to behave periodically in time. 
Thus one can take the solution in the form 

(3.10) 


where Z are magnitudes and A is the frequency (non-dimensional 
in this case since T 'is non-dimensional). 


Using expression (3,10), one can rewrite equations (3,9) 
as (dropping superscript (e) on Z for convenience) ; 

,a2_.2 2 A)[A©]+ (C 3 - C 3 , 2 X^- A^) 


[M] 


( Cl + C2 (A^( 1+Mc) - 202^. )§| 


dZ 


n= 0 



d^Z 

* 

- c., 


n= 0 

- ^^1 

d^Z 

drf 

+ ( ’X^( 1 +Mc ) 


d^Z 


C2 

mmrnmmmym 

drf ‘ 



> 


^ \ ct2 


Tl= H 


2 


1 ] = H 


(3,11) 
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Noting that t) = 0 and tj = H correspond respectively to 
nodes 1 and 2 of an element, one can rewrite equation (3.11) 
more briefly as: 



2 1 
^1 


1 

< 

PLKj+ Q [mb] + R [m]] 

Z’ 

^1 

= 

-Ml 

i 


^ < 

V2 ^ 




. M^ 


where 


P = 
Q = 
R = 



-C^ ((1 + Me) A. ^ - 2 QX) 

^3^^ ” C3. 2ft 

3 

“( + C 2 ( 1 Me) 


d^Z 

di1^ 



The above are the finite element matrix equations for 
a rotor element taking into account rotary inertia, shear 
deformation and gyroscopic effects. 


12) 



CHAPTER 4 


F.E. TRANSFER J/ATRICES AND METHOD OF SOLUTION 


1 Finite Element Transfer Atetrices 

The assembled matrix |^P[K] + Q [MB] + R [ M of 
equations (3.15) is a square 4x4 matrix. From equation 
(3,15) one sees that this assembled matrix relates nodal 
variables Z and Z* to the forces V and moments M at the 
nodes. To obtain the transfer matrix which provides a 
relation between variables at nodes 1 and 2, the equations 
(3.12) are rearranged. To achieve' this,the equation (3.1:^ 
are partitioned into submatrices as follows: 


'[A] 

[B]' 

< 

Wi' 

> = < 

> 

0 

1 M 

[D]_ 




(^2 


where 


^3 


and 



then one has 


[A] 


+ [ B ] \,W}2 = 

- IN 

[c] 

(wll 

+ [D] Wj = 

(B) 


(4.1 .1) 


1,2 , 


(4.1 .2) 
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From the first of equations (4.1.2) one gets 

= - Cb]-' [A] -[B]-’ 

Usingthis, from the second equation of (4.1,2) one gets 
1f} 2 - fc]- [DliB]-’ [A]] (.W}^ - ^ 

The above two equations are put together to yield the 
following transfer matrix equation: 

wl -iBr''i:A] -[B]-'' |[w}' 

[( 4)2 ' [CHDJCB]-' [A] -[D][B]-' |[F}J, 

( 4 . 1 . 3 ) 

Writing this more' concisely and for the general 
case of the jth element having nodes j and ( j'+l), 

= [T]j {,SV}^ (4.1.4) 

where (SV} . is the state vector at node j and [T] . is the 
field transfer matrix for the (finite) element j; 




31 


= [T]„ ... [T], Isv}, 

or 

^-^Wn+l = isv}, (4.2.1) 

where [OT] is the overall transfer matrix between state- 
vectors at nodes 1 and (n+1). Boundary conditions at 
these nodes can then be used to obtain a determinant of 
order two which goes to zero for natural frequency values 
of X • For critical speeds, forward whirl requires Q equal 
to X while backward whirl requires 0 equal to -x • 

4.3 Rotor with In-Span Supports 

In-span supports considered are all of the 'simple* 
kind allowing continuity of slope an-d bending moment while 
imposing the condition that deflection is zero, A clamped 
in-span support would be equivalent to having the shafts 
isolated from each other for the case of transverse vibrations. 

A simple in-span support introduces a shear disconti- 
nuity. It is assumed that other components of state 
vector {SV} remain continuous across such a support. It 
may be mentioned that this assumption implies an approxi- 
mation when considering shear effects on the shafts if a 
single governing equation is used because actually slope 
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due to bending moment alone would remain continuous. Since 
there is a shear disccantinuity, the slope due to shear effects 
and consequently the total slope will change across the support. 
However, since this effect is usually small it does not appear 
to be an unreasonable approximation. Difficulties with using 
a two equation approach where this distinction could be 
maintained have been pointed out earlier. 

Introduction of a shear discontinuity necessitates the 
derivation of a ’point transfer matrix’, briefly point matrix, 
for transferring state vector variables from one side of the 
support to the other. 

Two methods are available for handling such situations, 
Pestel and Leckie [9] follow a technique involving elimination 
of unknownsat different in-span supports and they reduce 
the size of matrices being handled by appropriate manipulations, 
Pilkey and Chang [10], on the other hand, derive point matrices 
of same order as field matrices using a condition from the 
station to the right of the in-span support. The latter is 
the more elegant method in terms of programming the computer. 
Subramanya [11] has used it in conjunction with transfer 
FEM, The same technique is used here. 

Two case s_ are handled: 

(a) in-span support with simple support on its right (to 

handle shafts with both extremities on simple supports), and 
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( b) in-span support with free end on its right (to take 
care of cases with overhang) 

In-span support with Simple Support on its Right 

The specific case of Figure 4.1(a) is considered with 
two in-span supports and simple supports at shaft 
extremeties. Generalising toa case of a n in-span supports 
will be a simple matter. 

This shaft has three spans, each simply- supported 
on both ends, field transfer matrix for span j 

while is the state vector at support point j; 

Since force V is discontinuous across in-span support, 
superscripts L and R are used on V to denote force just to 
the left and right respectively of support j. Other state 
variables being continuous across the support, no super- 
script is used on them. 

To obtain point matrix for support 2, transfer matrix 
relation for span 2 is used as follows: 


z 


"^11 ^12 "^13 ^14 


Z 

Z' 


• • • • 


Z’ 

1 



< 


M 

4 

P3 ! 

’^41 • * ^44 

, s2' * 

M 


(4.3.1) 
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Since = 0 one gets 


M2 7’ 

T.o S2 


where T. . are from T 

IJ S2 

Using the above expression, state vectors on the two 
sides of support 2 can be related as 


(4.3.2) 


tsv}«2 = [P]^ (SV}^2 

where [P] 2 is point matrix for support 2. 

Point matrix for support 3 can be obtained in a 
similar fashion using [T] 33 • overall relationship 


then becomes: 


'-SV}p4 [T]32i:Pl2 ms1 ^SV}, 


(4.3.3) 


or tsv}p4 - [OT] tsv}p, 

where [OT] is the overall transfer matrix, 



To get the determinant for finding the natural frequen- 
cies, along with unused boundary conditions of the system one 
is required to use the unutilized condition of an in-span 
support. For the present case, boundary condition = 0 
has been utilised in obtaining [P ]2 • From equation (4,3.3) 

using the remaining boundary conditions, one gets: 


OT 42 • * 0^43 • 'Vi “ ° (4.3.4a) 

while from transfer matrix relationsfor span 1 and making use 
of the unutilized condition Z ^2 ~ 0 ? 


12 ^p1 


13 ^p1 


= 0 


(4.3.4b) 


where "^^ 2 * "^13 from [t 3 • 

Putting equations (4.3.4) together one obtains the 
determinant to be satisfied as: 


”^12 ^13 


OT 42 0^42 


= 0 


(4.3.5) 


4.3.2 In-span support with Free End on its Right 

Consider the specific case of Figure 4, 1(b) . As before 
there are two in-span supports but point 4 is free unlike 
the previous case of section 4.3.1 . Point matrix for 
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Figure 4,1 : Multi-span rotors (a) both ends supported 

( b) with overhang 


support 2 is derived as in the previous section. For support 3, 
however, one needs to use a different condition, viz, = 0 

to derive the point matrix; 

Along lines similar to- those followed in the previous 
section, one has 

=[I],3tSV}«3 



Point matrix £F ]3 is hence: 




(4,3,6) 



Rest of the procedure is similar to that outlined in 


the previous section 4.3.1. 



CHAPTER 5 


RESULTS AND DISCUSSIONS 

FORTRAN programs have been developed for determination 
of natural frequencies and critical speeds using the one 
governing equation approach for a variety of cases. The two 
equation approach has also been used for the case of a uniform 
rotor . 

Results are presented in Table 1 through to Table 7 
and Figure 5.1 to Figure 5.2, The discussion that follov/s is 
with reference to particular tables and figures. Exact 
results where presented mean results obtained from a closed- 
form formulation. 

In the presentation of these results two values of factor 
Re (= k/L) have been used frequently. Small Re value (denoted 
by Re^ in some tables) of 0.0125 corresponds to a rotor with 
diameter being 1/20 of rotor length (’slender’ rotor) while 
large Re (denoted as Re^^) of 0,05 corresponds to a ’stubby' 
rotor with diameter being 1/5 of the total rotor length. For 
a steel shaft of circular cross section, material constant 
Me (= E/k'G) is 2.9 and this is the value used throughout. 
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Table 1 ; Critical speeds for simple-simple, cantilever and 
fixed-fixed rotors are shorn for various number of elements. 
Exact critical speeds available for simple-simple rotor [1] 
are also shown, , Study of this table shows that one can safely 
take 6 elements. For simple-simple case with 6 elements error 
in the first mode is 0.01?^ and in the fourth mode is 1,03?^. 

Some general trends that can be observed are: 

i) results converge from above, i.e, values obtained 
are larger than the exact values, and 
ii) , error is more in higher modes. 

These trends are usually observed in vibration problems 
solved using finite element method, 

Ficrures5.1 and 5,2 : These show the natural frequencies 

for different modes as functions of rotor speed, ^ * 

Figure 5,1 shows this for a uniform rotor for which shear 
deformation effects have been ignored. Natural frequency 
is seen to increase with rotor speed. Numerical difficulties 
did not allow determination of X for rotor speeds more than 
600 for the case considered (Re » 0.05). From a practical 
point of view.Q = 600 is a very high speed, being about 

sixty times the first critical speed. 

Figure 5.2 shows the natural frequency versus rotor 
speed plot for a simple-simple rotor with Re = 0.05 with shear 



effect included. Shear effect, it is seen, tends to reduce 
natural frequencies and also makes the curves flatten out 
at very high rotor speeds. 

In both cases, results are seen to match very well with 
those obtained from exact formulation [ • Hence the present 

technique can be used for quick and accurate determination of 
variation of natural frequencies for all types of rotor systems 
as a function of rotor speeds. 

Table 2 : Results in this case are from two governing 

equations approach. With exact integration, ill-conditioning 
was observed. Reduced integration was used to over come this 
problem [4,19 ] « Even with this, convergence is observed 
to be slow; For the simple-simple supports case, comparison 
with exact results shows that critical speed values converge 
to an answer that is less than the exact. 

Table 3 ; This gives results for three different rotors 
using one and two governing equations approaches for small 
and large values of Re, Problem of ill-conditioning of 
equations is seen to be serious for small Re values for which 
answers are seen to be' higher than those obtained from a 
single equation approach. For large Re, however, the opposite 


is seen to hold true 
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Due to the erractic behaviour of results due to ill- 
conditioning, this two governing equations approach has not 
been pursued further. 

A’ Results for three different cases of end supports 
are presented for both small and large Re values. For small 
Re, results are seen to be very close to the natural frequency 
values of an Euler beam while for large Re, these are signifi- 
cantly lower. This is to be expected since accounting for 
Timoshenko effects affects natural frequencies significantly 
for thick beams/ rotors. The effect of gyroscopic action of 

rotor elements is shown by a comparison of natural frequency/ 
critical speed values of Timoshenko beams and rotors. As 
seen from the table, gyroscopic effect tends to increase the 
natural frequency of the system and is quite significant for 
large Re and higher modes. This is to be expected since firstly 
large Re implies a thicker rotor, and secondly hi<^er modes 
imply larger curvature of shaft. Both imply large gyroscopic 
moment and hence increased effect of the gyroscopic action. 

Table 5 ; Some cases of multi-span beams/ rotors with 
different sets of end-supports have been considered here. 

The table also shows schematically the various cases considered 
With reference to these figures, the fractional numbers (just 
above the shaft line) give the non-dimensional span lengths whi 
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the whole numbers (just below the shaft line) give the number 
of elements used for the span. To check the validity of 
the method employed to handle such problems, results for Euler 
beams from the transfer FEM are compared with exact results [isl 
and are seen to be in good agreement. It may be mentioned that 
the ’exact* values quoted for the simple— simple case have been 
obtained from a graph in [18] while for the other cases these 
have been obtained from tables in [18] . Exact values [18] 
for the simple-simple case are therefore relatively inaccurate. 

Table 6 ; For a double-span beam, changing position of the 
in-span support changes the critical speeds of a rotor and 
this fact can have useful applications in systems where the 
designer has the freedom to select the position of this support. 
Table 6 gives the variation in critical speeds for some different 
positions of the in-span support for uniform rotors. 

Different sets of end-supports have been considered. These 
results v;ere seen to follow patterns similar to those obtained 
by Gorman [18] for Euler beams. 

For the uniform rotors considered, two observations 
can be made. Firstly, it is seen that critical speeds vary 
considerably as in-span support position is varied. Secondly, 
this variation is such that the difference between the first 
two critical speeds also varies. This difference is seen to be 
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maximum at about U = 0.3 for the symmetric cases (identical 
supports at both ends), at about ^ ^ 0.5 for the fixed-free 
case, and about }j. = 0.4 for simple-free case. The designer 
can therefore ensure that operational range of a rotor is well 
removed from critical speed regions by a judicious selection 
of in-span support position. 

A similar exercise can be carried out for other multi- 
span rotors. The advantage of having this efficient means of 
evaluating critical speeds for different configurations is 
obvious , 

Table 7 ; Results for two cases of non-uniform rotors (as 
shown in Table 7) are given. These are based on the character- 
istic value of Re^ . Values for rotors with small and large 
Re's are presented. As before, values for small Re are seen to 
be close to Euler beam results while there is significant 
change from these for large Re. 



CHAPTER 6 


CONCLUSIONS 

Transfer Finite Element Method has been applied 
successfully for determining natural frequencies and critical 
speeds of rotors in a variety of situations. This method offers 
the advantage of reduced computer memory requirement vis-a-vis 
standard Finite Element Method. It is hence- a useful technique 
for solving such problems in rotor dynamics. 
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TABLE 1 


Convergence Study for Various Rotors 
(One governing differential equation. Re = 0.05, Me = 2.9) 

1 . Simple- simple 



Number of elements 1 

Mode 

4 6 

8 

12 

Exact' L1 3 

1 

9.641 9.639 

9.638 

9.638 

9.638 

2 

36.014 35.905 

35 . 886 

35.878 

35.876 

3 

73.343 72.399 

72.224 

72.152 

72.135 

4 

122.897 113.415 

1 12.656 

112.342 

112.261 

2. Cantilever 



I 


: r 

1 

Number of elements 

Mode 

— 

4 

6 

8 

12 i 

1 

3.477 

3.477 

3.477 

3.475 

2 

20.460 

20.444 

20.441 

20.439 

3 

52.196 

51.951 

51 .902 

51.880 

4 

90.805 

90.088 

89.801 

89.676 


Continued. . . ; 
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TABLE 1 (Continued) ; 
3. Fixed-fixed 
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TABLE 2 

Convergence Study for Various Rotors 
(Two governing differential equations. Re = 0.05, Me = 2.9) 


1 . Simple-simple 


1 



Number of 

elements 



6 

9 

12 

15 

18 

Exact [ 


9.746 

9.570 

9.510 

9.482 

9.467 

9.638 

■ 

37.960 

35.569 

34.779 

34.423 

34.231 

35.876 


83.398 

73.479 

70.307 

68.895 

68.142 

72.135 

■ 

144.764 

120.582 

112.599 

109.072 

107.201 

112.261 : 


2. Fixed-fixed 





Number of 

elements 

' 1 

Mode 

6 

9 

12 

15 

18 

1 

20.334 

19.586 

19.335 

19.221 

19.159 

B 

52.381 

48.126 

46.735 

46.109 

45.774 

B 

97.434 

85.319 

81 .354 

79.583 

78.635 

■ 

153.792 

129.629 

121.149 

117.348 

115.319 









Results from Transfer FEM 



10 Q 100 1000 10000 

FIG 5.1 Natural Frequency A vs Rotor Speed n 
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TABLE 3 

Critical Speeds from one and two Governing Equations Approach 
(Re^ = 0.05, Reg = 0.0125, Me = 2.9) 

Notes: i) Columns numbered 1 and 2 contain results from one 

and two governing equations respectively 

ii) 8 elements considered for one equation approach 
while 12 eleiiBnts considered for two equation 
approach 


1 . Simple-simple 



ReL 

— 

Reg 

Mode 

1 

2 

1 

2 

1 

9.638 

9.510 

9.855 

9.925 

2 

35.886 

34.779 

39.257 

40.336 

3 

72.224 

70.307 

87.767 

93.374 

4 

112.656 

112.599 

154.811 

172.535 


2. Cantilever 



, 

RSl 

Reg 

Mode 

1 

2 

- - - •- 

1 

2 

1 

3.477 

3.442 

3.513 

3.514 

2 

20.441 

19.508 

21.931 

22.191 

3 

51.902 

48.792 

61 .030 

63.244 

4 

89.801 

85.119 

118.620 

127.508 


Continued, 
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TABLE 3 (Continued) : 

3. Fixed-fixed 



: 

■ 

Reg 

Mode 

1 

2 

1 

2 

1 

21 .660 

19-335 

23.336 

22.580 

2 

54.682 

46.735 

61 .286 

63.520 

i 

3 ■ 1 

94.827 

81 .354 

119.392 

128.381 

4 

136.970 

121 .149 

195.836 

221 .439 
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TABLE 4 

Comparison of Critical Speeds for Small and Large Re Values 
(Reg = 0.0125, Rej^ =0.05, Me = 2.9, 8 elements) 


1 . Simple-simple 


Mode 

Euler 

[18] 

Reg 

ROl 


Rotor 

Timoshenko 

Rotor 

1 

9.87 

9.84 

9.85 

9.43 

9;64 

m 

39.48 

39.02 

39.25 

33.81 

35.89 


88.83 

86.65 

87.76 

66.56 

72.22 

1 

157.90 

151 .53 

154.81 

103.51 

112.66 


• 2. Cantilever 


lode 

Euler 

[181 

Re^ 

3 

Rej 


Timeshenko 

Rotor 

1 

Timoshenko 

Rotor 

1 

3.52 

3.51 

3.51 

3.44 

3.48 

2 

22.03 

21.82 

21.93 

19.25 

20.44 

3 

6147Q 

60.34 

61 .03 

47.33 

51.90 

4 ' 

120.91 

116.25 

118.62 

81.09 

89.80 


Continued, . ; . ; 
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TABLE 4 (Continued) : 

3. Fixed-fixed 


3de 

Euler 

[18] 

Re, 

3 

} "■ 

! 

Re 

L 

Timeshenko 

Rotor 

t 

Timoshenko 

Rotor 

1 

22,21 

22.29 

23.34 

21 .22 

21.66 

2 

61 .67 

60.87 

61 .29 

52.11 

54.68 

3 

120.91 

1 17.75 

119.39 

89.90 

94.83 

4 

199.86 

191 .55 

195.84 

131.08 

136.97 
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TABLE 5 

Natural Frequencies/Critical Speeds for Some Multi-span Cases 

(Re = 0.05, Me = 2.9) 

0.3 0.7 

Double-soan: 

r- ye 71; A 

2 3 


Mode 

j 



Euler beam 

Rotor 

Exact [is] 

Present 

1 

26.01 

26.37 

24.95 

2 

85.56 

87.09 

70:62 

3 

138.06 

140.75 

104.22 

4 

201 .64 

228.52 

138.36 


Triple-span: 


0.4 0.3 0.3 


a: 


X 


X A 


4 2 2 


Mode 

Euler beam 

Rotor 

Exact [18] 

Present 

1 

26.23 

26.23 

23.88 

2 

107-y54 

107.60 

86.40 

3 

161.80 

162.76 

111.91 

4 

228.61 

230.95 

142.26 

Continued., 
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TABLE 6 

Variation in Critical Speeds of Double Span Rotors with 
Change in Position of In-span Support 


yv TT 



(0 

11 

0.0125) 





p- 

Mode 

0.1 

0.2 

0.3 

0.4 

0.5 

1 

17.83 

21 .28 

26.25 

33.29 

39.25 

2 

57.80 

69.75 

84.99 

76.09 

61.24 








^ ^ 




Mode 

0.2 

0.4 

0.5 

0.6 

0.8 j 

1 

1 

4.90 

7.53 

9.85 

13.52 

21 .83 1 

2 

31.13 

49.94 

60.95 

51,87 

51.82 i 


Continued 
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TABLE 6 (Continued) : 


^ ^ ; 15, 




Mode 

0.1 

0.2 

0.3 

0.4 

0.5 

1 

26.91 

31 .91 

39.53 

50.86 

61.24 : 

2 

72.20 

88.09 

109.24 

111 .18 

88;80 : 


A Ts: 

h— ^ 



: 

Mode 

0.2 

0.4 

0.5 

0.6 

0*8 1 

1 

4.75 

7.05 

9.05 

11.98 

14.57 

2 

30.46 

47.00 

46.24 

36.65 

44.96 









TABLE 7 

Critical Speeds for Non-uniform Rotors 
(Results based on characteristic value of Re^) 


I* — • 0;5 


Re 2 Re^ R.e^ = 0.025, Re 

1 

— ^ 

*\ Rot^- Re. = 0.05, Re 


Mode 


Euler 


10.89 

44.35 

98.32 


10.80 
43.05 
91 .90 


10.56 

39.35 

76.31 


0.25 


Re^::^ Re^=Re2=0.025,Re2=0.C 
Rej^=^ Re^=Re2=0,05, Re 2=0. ( 


Mode 

Euler 

Reg 


1 

11.75 

11.67 

11.44 

2 

43.55 

42.26 

38;62 

3 

96; 85 

90.31 

74; 77 
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APPENDIX A 

NON-DIME NSIONALISED GOVERNING DIFFERENTIAL EQUATIONS 


Characteristic length : L 

Characteristic time, T* = ( mL'^/EI) where m = p A 
Non-dimensional variables used are therefore: 

Z = z/L, S = s/L, T = t/T*,Q = to 

Then 


3^z ^ 

1 

3^2, 


3^z ^ 

L 

3 s^ 

L(n,1) • 

3 

♦ 

St"^ 

_n * r~Zn 


1 3? 




1 3^ 

"ds 

t • rs 


f 

St”^ 

T* *3 1’^ 


A1 . Two-qoverninq Differential Equations 
Equations (2,2.3) are; 
mSr - ^ [k'GA ( z* - y )] = 0 


3s 


(El + k’GA ( z’ -W ) -mk 


2 


+ 2imk (i) ? = 0 


Using non-dimensional variables and expressions the 
above equations become: 


- S M - i ^ [>='QA(|| -n] = 0 


3Z 


t; 3 T 

1 , 8 , El 


- T • ^ (5i - ft- ^ 


^ 3% . mk^n 3^ 


+ 2i Q 


y2 ■ 3T 



A 


For a uniform rotor section, k'GA and El are constants. 
The above equation can hence be rewritten as: 

(z" - W ' ) - Mc.Re^ . Z = 0 

^ " + 1 — ^ - Re^ (i* - 2iQ '? ) = 0 

Me. Re 

where Me is the non-dimensional material constant, E/k’G, 
and Re (also non-dimensional) is the rotary inertia factor, 
k/L, Me accounts for shear effects while Re accounts for 
rotary inertia of elements and also contributes to the 
gyroscopic action term (since polar moment of inertia is taken 
to be twice equatorial moment of inertia for a circular disc 
for equation (2.2.3)). 

The above is presented as equation (2,2.4) in Chapter 2, 

A, 2. Single- governing Differential Equation 
Equation (2,3,3) is 


El 

m 


.IV 


k^ ( 1 



• • !I 

Z 


+ 2i o k^ 


•n 

Z 


+ Z + 


mk^ •; 
k'GA ^ 


9 

•2io>A z = 0 


Using non-dimensionalised parameters one obtains: 



A. 3 


El 

d^Z 

- ( 1 + 

k'G^ * 

Cl** 

1 3 Z 



mL'^ 

• “Z 

as 

ltJ as^T^ 



2i Q 



L 

3^z mk^ 

L 

a^z 

k — ^ 

- ltJ 

2 

as ai 

j2 

^ * 

* k'GA • 

4 

a 


- 2i .2 

mk^ L 

3 3^ 

= 0 . 




k'GA ^4 

^ * 

ai'^ 




2 s 

Dividing the equation by ( L/T and simplifying 
expressions: 

Z ^ - Re (1+Mc) Z + 2i.a. Re . Z + Me. Re . Z 

+ Z - 2i Q , Ate.Re"^ . *Z* = 0 

with Me and Re as defined above (See. A.1) . This is presented 
as equation (2.3.4) in Chapter 2. 



APPENDIX B 


3,1 


FINITE ELEfvENT FORMULATION FOR TWO GOVERNING EQUATION APPROACH 


For a finite element, governing equations from equation 
(2.2.4) are: 

(Z" -? •) - . 2 = 0 (B.1a) 

's" + (Z* -^) - (V - 2i0¥) = 0 (B.lb) 

where 

’ : derivative w.r.t, local non-dimensional coordinate 

(and not w.r.t. S as in equation (2.2.4)) 

; ; derivate w. r. t. non-dimensional time T 

= Mc.Re^, = Re^ 

Z = 2 (■/!], T),? = ^' (R, T) . 

Assume polynomial solutions for Z and ? over the eleiaent: 

Substituting the above expressions in equations (B..1) 

Will give residues from both equations. Minimising residue 
from (B.la) by Galerldn’s method using iN^J and that from (B. 1b’ 
using |_N J, one gets: 



H 




H 


- C^ ] {n^} LN"i dTl(2'}^““^ = tN^}l 


. - ( ne) 


.H 


z\x« (e) 


- {N^} 


32(e)l H 


3t1 


(B.2a 


H 


I / 


H 




1 b 


J_ 

C. 


H 


I [N'^}\N"jdn 


^ - H" 1 I . ."V| . - • • ^ ( ^ ^ . 


“ ^2 f ^N ^ \_N 1 dTl (w ] + 2i C 2 J } [N J dt) f 


= - 


3f (e) 


3T1 


H 


(B.2b 


Applying the criteria of completeness and compatibility 


it is seen that interpolating functions for both Z and W need to 

rZ- 


linear. In terms of shape functions, [N ] and [ N ] therefore 
become identical and these are [12] : 

LnJ = t 1 -?i § J 

ne) 


while 


12, 


=L5^ y 


and 



where subscripts 1 and 2 denote values at nodes 1 and 2 of the 
element (i.e. = 0 and H respectively). 

Assuming and of e^^”^ form, one gets 

from equations (B.2): 

[K.1][Z]^’^®^ = (N} ^ 

- (n) ^ 

o ( B •Sc 

[K3] ^ [ K5] [K4]£W}^’^®^ 

-C 2 . [f«] + 2 Q.X ; [MM]^^ = N 

(B. 3 : 

where 

H , , H , 

[K1] = ]■ .[N] iNldtl , [K2] = ]• (.N}\_NJ d^ 

o o 

H - H 

[K3],= J (n'} [^n'] dn = [K1] , [K43 = XtNjLNjdTl I 

o o 

H H ^ 

[ K5] = S (N} \.n’] d^ , [MM] = X W iNj dt) i 

O '.O . 

are all square matrices of order 2. Note that while [ Kl] and ' 
are identical, a distinction has been maintained between [K4]4 
[ MM ], This because due to ill-conditioning, reduced integral 
by Gaussian quadrature is carried out for[K4] , The elements 
of these are given in Appendix C, 



Equations (B.3) can be written compactly as: 


[[K1 ]- Ci . - [K2j{^}^ 


t'oi f'P ■)( -- 


where 


(B.4; 


- -^ [ K5] + |^[K3] + [K4] + (2CX_'x2^ j. 


( ne) 


(B.4 


-?) , and M. = 


d li 


The above equations can be combined to yield the followin 
finite element matrix equation; 


[K1] - X [KW] 


-[K2 3 


- [K5] [K3] +[K4] [MM] 



APPENDIX C 


C.1 


Matrices 


[K] 


H' 


£M] = 


H 

420 


[MB]= ^ 


FINITE ELEMENT ?4ATRICES 

for four degrees of freedom finite elements 


12 

6H 

-12 

6H ■ 


CM 

-6H 

2 

2H"^ 

symmetric 


12 

-6H 

- 



4H^ 

156 

22H 

54 

-13H 


2 

4H 

13H 

2 

-3H^ 

symmetric 


156 

-22H 

- 



4h2 

6 

H 

6 

H 

5 

10 

5 

10 


2 

2H^ 

H 



15 

10 

" 30 

symmetric 


6 

5 

H 

10 


2 

2H^ 

15 


C.2 




v-' a 


Matrices for two degree of freedom element: 


[Kl] = [K3] = ^ 

[ K2] = -i 

[K4]= I 
[K4]= § 

[ K5 ]= ^ 

[ MM] = I 


1 

-1 

-1 

1 

2 

1 

1 

1 

-1 



2 

1 


-1 


1 

-1 

1 

1 

2 

1 

1 

l' 

1 

r 

2 


(without reduced integra 


(with reduced integratio 
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